package com.xy6.algo.divide;

/**
 * 矩阵乘法
 * 
 * <pre>
 * 请编程实现矩阵乘法，并考虑当矩阵规模较大时的优化方法
 * 两个矩阵的乘法仅当第一个矩阵B的列数和另一个矩阵A的行数相等时才能定义。
 * 如A是m*n矩阵和B是n*p矩阵，他们的乘积C=AB是一个m*p矩阵，它的一个元素（i,j）满足
 * 1=<i<=m, 1=<j<=p
 * 计算公式：
 * C[i][j] = sum(A[i][k] * B[k][j])  k属于[1,n]
 * 时间复杂度: m*n*p
 * 如果A、B维度都为n，则时间复杂度为O(n^3)
 * 
 * Strassen算法将相乘次数由8次降为7次，时间复杂度为O(n^2.8)
 * </pre>
 * 
 * @author zhang
 * @since 2018-01-21
 */
public class MatrixMultiply {

	public static void main(String[] args) {
		int[][] A2 = { { 1, 2 }, { 3, 4 } };
		int[][] B2 = { { 0, 1 }, { 0, 0 } };
		int[][] C2 = multiply(A2, B2);
		print(C2);
	}
	
	/**
	 * 通用矩阵相乘算法
	 * 
	 * <pre>
	 * 时间复杂度: m*n*p
	 * </pre>
	 * 
	 * @param A
	 * @param B
	 */
	public static int[][] multiplyForce(int[][] A, int[][] B) {
		int m = A.length;
		int n = A[0].length;
		int p = B[0].length;
		int[][] C = new int[m][p];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < p; j++) {
				for (int k = 0; k < n; k++) {
					C[i][j] += A[i][k] * B[k][j];
				}
			}
		}
		return C;
	}

	/**
	 * 打印一个矩阵
	 * 
	 * @param C
	 */
	public static void print(int[][] C) {
		for (int i = 0; i < C.length; i++) {
			for (int j = 0; j < C[0].length; j++) {
				System.out.print(String.format("%d ", C[i][j]));
			}
			System.out.println();
		}
		System.out.println("-------------");
	}

	/**
	 * Strassen算法求解矩阵相乘
	 * 
	 * <pre>
	 * 适应范围：矩阵A、B行数、列数相同，且为2的整数倍
	 * 时间复杂度: O(n^lg(7))
	 * http://www.mamicode.com/info-detail-673908.html
	 * 蛮力法求解进行8次乘法，Strassen算法求解进行7次乘法
	 * O(n^3) vs O(n^lg(7))，随着n的变大，比如当n >> 100时，
	 * Strassen算法是比通用矩阵相乘算法变得更有效率。
	 * </pre>
	 * 
	 * @param A
	 * @param B
	 */
	public static int[][] multiply(int[][] A, int[][] B) {
		if (A.length <= 1 || B.length <= 1) {
			int[][] C = new int[A.length][B.length];
			C[0][0] = A[0][0] * B[0][0]; 
			return C;
		}

		int n = A.length;
		int[][] a = subMatrix(A, 0, 0, n / 2 - 1, n / 2 - 1);
		int[][] b = subMatrix(A, 0, n / 2, n / 2 - 1, n - 1);
		int[][] c = subMatrix(A, n / 2, 0, n - 1, n / 2 - 1);
		int[][] d = subMatrix(A, n / 2, n / 2, n - 1, n - 1);

		int[][] e = subMatrix(B, 0, 0, n / 2 - 1, n / 2 - 1);
		int[][] f = subMatrix(B, 0, n / 2, n / 2 - 1, n - 1);
		int[][] g = subMatrix(B, n / 2, 0, n - 1, n / 2 - 1);
		int[][] h = subMatrix(B, n / 2, n / 2, n - 1, n - 1);

		int[][] p1 = multiply(a, minus(f, h));
		int[][] p2 = multiply(plus(a, b), h);
		int[][] p3 = multiply(plus(c, d), e);
		int[][] p4 = multiply(d, minus(g, e));
		int[][] p5 = multiply(plus(a, d), plus(e, h));
		int[][] p6 = multiply(minus(b, d), plus(g, h));
		int[][] p7 = multiply(minus(a, c), plus(e, f));
		
		int[][] t1 = plus(minus(plus(p5,p4), p2),p6);
		int[][] t2 = plus(p1,p2);
		int[][] t3 = plus(p3,p4);
		int[][] t4 = minus(minus(plus(p1, p5), p3), p7);
		
		return join(t1,t2,t3,t4);
	}

	/**
	 * 将四个子矩阵合并为一个矩阵
	 * 
	 * @param t1
	 * @param t2
	 * @param t3
	 * @param t4
	 * @return
	 */
	private static int[][] join(int[][] t1, int[][] t2, int[][] t3, int[][] t4){
		int[][] t = new int[t1.length + t3.length][t1[0].length + t2[0].length];
		int row1 = t1.length;
		int col1 = t1[0].length;
		for(int i=0; i<row1; i++){
			for(int j=0; j<col1; j++){
				t[i][j] = t1[i][j];
			}
		}
		for(int i=0; i<t2.length; i++){
			for(int j=0; j<t2[0].length; j++){
				t[i][col1 + j] = t2[i][j];
			}
		}
		for(int i=0; i<t3.length; i++){
			for(int j=0; j<t3[0].length; j++){
				t[row1 + i][j] = t3[i][j];
			}
		}
		for(int i=0; i<t4.length; i++){
			for(int j=0; j<t4[0].length; j++){
				t[row1 + i][col1 + j] = t4[i][j];
			}
		}
		
		return t;
	}
	
	/**
	 * 求一个矩阵的子矩阵
	 * 
	 * @param A
	 * @param x1 起点行号
	 * @param y1 起点列号
	 * @param x2 终点行号
	 * @param y2 终点列号
	 * @return
	 */
	private static int[][] subMatrix(int[][] A, int r1, int c1, int r2, int c2) {
		int[][] a = new int[r2 - r1 + 1][c2 - c1 + 1];
		for (int i = r1; i <= r2; i++) {
			for (int j = c1; j <= c2; j++) {
				a[i - r1][j - c1] = A[i][j];
			}
		}
		return a;
	}

	/**
	 * 两个矩阵相加
	 * <pre>对应元素相加</pre>
	 * 
	 * @param A
	 * @param B
	 * @return
	 */
	private static int[][] plus(int[][] A, int[][] B) {
		int[][] C = new int[A.length][A[0].length];
		for (int i = 0; i < C.length; i++) {
			for (int j = 0; j < C[0].length; j++) {
				C[i][j] = A[i][j] + B[i][j];
			}
		}
		return C;
	}

	/**
	 * 两个矩阵相减
	 * <pre>对应元素相减</pre>
	 * 
	 * @param A
	 * @param B
	 * @return
	 */
	private static int[][] minus(int[][] A, int[][] B) {
		int[][] C = new int[A.length][A[0].length];
		for (int i = 0; i < C.length; i++) {
			for (int j = 0; j < C[0].length; j++) {
				C[i][j] = A[i][j] - B[i][j];
			}
		}
		return C;
	}

}
